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We investigate the early part of the crystal nucleation process in the hard sphere fluid using 
data produced by computer simulation. We find that hexagonal order manifests continuously in 
the overcompressed liquid, beginning approximately one diffusion time before the appearance of the 
first ‘solid-like’ particle of the nucleating cluster, and that a collective influx of particles towards 
the nucleation site occurs simultaneously to the ordering process: the density increases leading to 
nucleation are generated by the same individual particle displacements as the increases in order. 

We rule out the presence of qualitative differences in the early nucleation process between medium 
and low over compressions, and also provide evidence against any separation of translational and 


orientational order on the relevant lengthscales. 

I. INTRODUCTION 

As a typical first-order phase transition, crystallization 
from the metastable melt begins with a nucleation pro¬ 
cess. There are two order parameters which characterize 
the transition (density and ‘crystalline order’). It is not 
a priori evident that both parameters undergo the same 
dynamics during the transition process. Consequently, 
a “density-first”, a “bond-order-first”, and several more 
complex phase transition scenarios have been proposed 
and vividly discussed in the literature. 

In this article we investigate the role of density and 
order fluctuations prior to nucleation and in the envi¬ 
ronment of the growing nucleus. As a model system we 
use spherical particles that have repulsive interactions 
only. A hard sphere model is the starting point for 
many theoretical treatments of granular, fluid, glassy 
and crystalline systems, and may be sufficient without 
further refinement if excluded volume interactions are 
more significant than long-range forces. Technological 
examples where a hard sphere model is sufficient within 
important regimes include the study of metal solidifi¬ 
cation (e.g. [1]) and the formation of colloidal crystals 
(e.g. [2]). Crystallization in hard spheres has been 
studied extensively but retains many open questions, 
of which the nature of the initiating fluctuation is a 
particularly active concern. 

Classical Nucleation Theory (CNT) assumes that the 
emerging nucleus and the surrounding fluid possess the 
properties of the respective bulk phases, i.e. that the 
emerging nucleus already has the order and density of 
the bulk crystal and that the surrounding fluid is not 
affected by the growth of the nucleus. This picture is not 
immediately credible in real liquids and colloidal suspen¬ 
sions, where transport of material is relevant. A detailed 
consideration of the outcome of mass conservation and 
finite transport speeds in the fluid has been made for 
vapour-liquid nucleation by Lutsko et al. [3-5], leading 
to an expectation of substantial early density changes 


on a lengthscale larger than the initial nucleus (for 
the vapor-liquid nucleation of Lennard-Jones particles). 
Early densification in fluid-solid transitions is not a 
priori expected to be as important as in vapour-liquid 
transitions, due to the smaller difference in densities 
of the two phases, however this phenomenon has been 
reported in computer simulations of hard sphere crys¬ 
tallisation by Schilling et al. in 2010 [6] as well as from 
experiment [7, 8]. 

A computational study by Russo & Tanaka [9] has 
examined structure and density changes for a set of 
nucleating trajectories at a number density ~ 1.02, 
where a is the diameter of a particle. Augmenting 
these trajectories with fluctuation data drawn from 
the metastable liquid, they note a coupling of order 
and density, but state that “the density increase is 
foreshadowed by the prestructuring of the nucleus”, 
adopting a position which we will crudely summarise as 
‘order-first’. 

Tan et al. presented, from optical microscopy, two 
distinct nucleation pathways for their colloidal system. 
Both pathways began with hexagonally-ordered precur¬ 
sors then developing into either bcc or fee metastable 
structures. A three step process was described, with 
the final phase being either bulk bcc or rhep depending 
on the radius of the particles. The presence of bulk 
bcc, and also the dependence on radius, indicate that 
forces other than excluded volume were important in 
these experiments, however they remain an interesting 
reference in that the authors state that the nucleation 
sites were not correlated with increased local density. 
In fact the authors claim that “nucleation rarely starts 
from the densest regions” [10], see also a comment 
on this work by Granasy & Toth [11]. Another in¬ 
vestigation using microscopy, by Lu et al. [12], also 
found a decoupling of the density and the nucleation 
event, stating: “nucleation events were observed that 
rarely start from the denser regions of colloidal samples”. 
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Kawasaki & Tanaka [13] found medium range ordered 
precursors of sizes larger than the critical nucleus in 
simulations of hard spheres. The lifetime of these 
structures was estimated as “a few times the relaxation 
time Tc”. It was stated that the ordering was not icosa- 
hedral but rather hep with multiple defects. A similar 
observation was reported by Schilling et al. [6, 14] where 
low symmetry clusters of considerable size spontaneously 
transformed into highly ordered crystals. The effect of 
ordered precursor structures on crystallization should be 
to reduce the surface tension, which would be consistent 
with experiments of Gasser et al. [15] finding strongly 
aspherical shapes for the critical or near critical nuclei. 
The most recent minireview [11] states that fluid-solid 
nucleation is probably multistep in most cases, but that 
the specific steps expected to be relevant in a given 
system are to date unclear. 

The overview of the current literature is confusing: 
is the now-common idea of hard sphere nucleation as 
a two step process physically useful? Are densification 
and ordering separable from each other, and if so does 
one in general come first? The emergence of unclear and 
apparently contradictory statements creates a need to 
re-examine the onset of crystallization. It is not desir¬ 
able to reduce nuanced expositions from the literature 
to blunt and simplified statements, however we must 
remark that the ‘density first’ [6] and ‘order first’ [9] 
positions cannot both be right. The summarised claim 
in Schilling et al. 2010 [6] is that of a two step mech¬ 
anism, remarkable in that it arises without attractive 
forces: “The metastable fluid relaxes the density first, 
by producing dense low symmetry clusters”. The key 
claim from Russo & Tanaka 2012 [9]) is of the opposite 
two step mechanism “The transition from liquid-like to 
crystal-like happens at constant density”. Kawasaki & 
Tanaka argue in a separate work [13] that the sequence 
in detail is liquid hep dense rhep. 

To clarify the sequence of events we argue that, at 
least at low and medium over compressions, the fluid 
does indeed possess correlations over medium length- 
scales in hexagonal order however these are typically 
accompanied by synchronised density fluctuations of 
the same lengthscale: the initiating fluctuation is of 
medium range and of low amplitude in both density 
and order, and has the same pattern of radial decay in 
each parameter as for the normal quiescent fluctuations. 
Our analysis shows a synchronised increase of density 
and hexagonal order leading up to nucleation, although 
we allow that elevated hexagonal order is probably the 
more useful (less noisy) of the two reaction coordinates 
as a predictor of nucleation events. 

A further issue which has been mentioned in relation 
to this subject is the existence of two distinct ordered 
phases in 2D systems of hard discs. The hard disc phase 
intermediate between the fluid and solid (called the 


hexatic) is distinct from the crystal in that is has long- 
range orientational order and short-range translational 
order [16]. There is no known hexatic-like phase of 
bulk 3-spheres, however the possibility of some related 
unstable state manifesting along the nucleation pathway 
has generated a certain amount of excitement, therefore 
translational and orientational order are typically dis¬ 
tinguished in the literature in case they should turn out 
to be meaningfully different to each other. 

We argue that to claim a separation of translational 
and orientational ordering in 3 dimensional hard spheres 
would require much stronger evidence than has yet been 
seen: well-defined and distinct scaling of orientational 
and translational order in the fluctuations prior to 
nucleation events would need to be shown, and this 
has not been observed. Orientational order is conven¬ 
tionally defined with respect to the neighbourhood of a 
single particle via a decomposition of the orientational 
distribution of ‘bonds’ (displacement vectors between 
neighbouring particles) into spherical harmonics {qe^Ve 
etc.), while translational order is typically computed us¬ 
ing reciprocal-space observables not directly comparable 
to the orientational analysis. We attempt to test this 
distinction here by comparing translational and orien¬ 
tational order as directly as possible, via a direct-space 
analysis of bond lengths which is more comparable to 
the spherical harmonic analysis of bond orientations 
conventionally used to report the orientational order. 

In this work we analyze new simulations at low volume 
fractions and using different dynamical schemes, and we 
also revisit the 2010 dataset. We show that local and 
also medium-range density changes take place simulta¬ 
neously to the initial formation of (weak) translational 
and orientational order, presenting an order-with-density 
model which is distinct from the order-first mechanism 
preferred by Russo & Tanaka [9]. We also test the 
expectation based on work by Kawasaki et al. [13], 
Tan et al. [10], Barros & Klein [17], and Schilling et 
al. [6, 14] that weakly ordered precursors to crystalliza¬ 
tion should or might be present, and And that hexagonal 
ordering manifests gradually and continuously in the 
fluid prior to the formation of the first definitively solid- 
like particles, without evidence of any intermediate state. 


II. RE-ANALYSIS AND EXTENSION OF THE 
2010 NUCLEATION DATASET 

As outlined in the introduction, previous studies of 
simulation data have suggested an important role for 
dense precursors in the structure formation process of 
hard spheres [6]. We revisit an extant dataset com¬ 
prising independent runs of N = 216000 hard spheres 
of diameter a at a number densities of = 1.03, 
1.027 and 1.024, corresponding to chemical potential 
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Scheme 

N runs 

N particles 

[ksT] 

Year 

1.030 

MC 

4 

216000 

-0.585 

2010 

1.030 

MD 

6 

216000 

-0.585 

new 

1.029 

flatPERM-FFS 

21 

20262 

-0.58 

new 

1.027 

MC 

6 

64000 

-0.56 

2010 

1.012 

flatPERM-FFS 

15 

19924 

-0.50 

new 


TABLE L Simulations discussed in this work. Runs using 
rare event sampling (flatPERM-FFS) were computationally 
cheaper than brute force MD and MC, therefore more runs 
were made (and smaller errorbars achieved). 


differences in of about |A/r| O.G/c^T per particle. 
Nucleation at these over compressions was observable 
without the use of accelerated sampling algorithms. 
The time evolution was defined by MC simulation 
mimicking Brownian diffusion, with particles undergoing 
independent displacements at a fixed attempt frequency 
and small random stepsize much less than a. In order 
to extend the original datasets, further runs were 
made, using event-driven molecular dynamics (MD) at 
= 1.03 thus testing for effects arising from the choice 
of dynamical scheme, and also using MD simulations 
accelerated by fiat histogram Pruned and Enriched 
Rosenbluth Forward Flux Sampling (fiatPERM-FFS) at 
lower overcompressions (tab. I). 


The known sampling difficulties for FFS calculations 
in this system [19] were ameliorated using the fiat 
histogram Pruned Enriched Rosenbluth FFS method 
(fiatPERM-FFS) [20-22] which was added to the func¬ 
tionality of the freshs .org [27] sampling system for the 
purpose of this calculation. flatPERM-FFS was found to 
give significant improvement in convergence relative to 
direct FFS. flatPERM-FFS is a path sampling technique 
applicable to stochastic dynamical systems subject to 
rare events, i.e. bottlenecks in their dynamics such as 
crystal nucleation. The scheme operates by selectively 
branching multiple copies of trajectories which make 
progress with respect to a specific collective variable 
(here the number of particles in the largest crystalline 
cluster). As the dynamics are stochastic, the branched 
trajectories diverge with a proportion of them making 
further forward progress in the reaction coordinate. 
The gist of the gain from using fiatPERM rather than 
direct FFS is that paths are selected for branching in a 
Bayesian way that takes account of their histories. 


Voronoi volumes were determined using Rycroft’s 
voro++ tool [18]. The crystalline cluster was identified 
by the dot product of Steinhardt bond order parameter 
qe of neighbouring particles [23-25]: 


= ^ E (1) 

qeqeihj) = ^ E C(^)C*(j) (2) 

m =— 6..6 

where Y^(f) is the of 13 complex components of 
the sixfold spherical harmonic, and where a bond was 
treated as solid-like for the purpose of analysis if the 
ratio qeqeii, j)/^/q6q6{i,i)q6q6{j, j) was greater than 0.7. 

Hexagonal/icosahedral order was also measured via the 
third-order parameter wq: 

= E (ri n o) (3) 

m,n,o V / 

s.t. m-\-n-\-o=0 

where parentheses indicate a Wigner 3-j coefficient 
determined for the m^n^o. wA^qA and wQ bond order 
parameters were found using the bop utility [28]. In 
all plots, the g6 value is calculated over the 12 nearest 
neighbours j for each particle, in order to have the most 
sensitive probe of hexagonal order. 

The progress coordinate of the flatPERM-FFS calcu¬ 
lation was chosen as the number of solid-like bonds (over 
a threshold of 0.75) in the largest cluster of particles 
having 10 or more such bonds. The initial interface A 
was defined as 20 bonds, with further interfaces placed 
not less than 20 bonds apart. The interface placement 
algorithm of Kratzer et al. [26] was used. Only one 
trajectory per FFS calculation was analysed, in order to 
avoid complicating the calculation of errorbars by the 
use of partly correlated trajectories. The University of 
Luxembourg HPC service was used [29]. 

In order to analyse the local environment of particles 
from an additional perspective, the parameter-free algo¬ 
rithm of van Meel et al. [31] was used to estimate the 
number of actual neighbours per particle, with a value of 
12 strongly indicating either icosahedral or close-packed 
order. This algorithm operates by associating a solid 
angle to each additional neighbour moving outwards 
from i until the solid angle subtended is Air. 

In fig. 1(a) we show the average size of the largest 
crystalline cluster. The definition of solid particles for 
the purpose of plotting this figure was chosen to be the 
same as in ref. [6], such that particle i was defined as 
solid if 11 or more neighbours j gave a g6g6 product 
> 0.7. This differs slightly from the thresholds chosen 
to define cluster membership in order to construct a 
reaction coordinate for the FFS runs, of 10 and 0.75. 
Insensitivity of this statistic to threshold choice is shown 
in the supplementary data (fig. SI [30]) by presenting 
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FIG. 1. (a) Number of particles in the largest crystalline clus¬ 
ter present, (b) g6 hexagonal order statistics (c) Voronoi vol¬ 
ume, (d) per-particle translational order given as the standard 
deviation of neighbour distances within the first and second 
shells, and (e) the number of neighbours found using the van 
Meel algorithm. Averages are over a sphere of radius 2cr (see 
text). Dashed lines show the whole-system average Voronoi 
volume. Errorbars are twice the estimated standard error. 


the same statistic for a minimum cluster size of 7, as in 
[9]). Because nucleation events are randomly distributed 
in time, it was necessary to define a different time-zero 
for each trajectory such that the time series collected 
over different nucleation events could be averaged 
together. A space-time zero marking the centre and 
start of each nucleation event was defined. The position 
of the nucleation event was defined as the centre of 
mass of the largest cluster in the system at the final 
timepoint when this cluster was of size 30 particles 
(i.e. at the time after which it grew irreversibly). The 
time-zero for the nucleation event was then defined by 
following the 30 particles backwards in time until all of 
them were liquid-like. The time zero is thus defined as 
the first final solidification of any particle in this clus¬ 
ter. Averaged quantities were then taken over particles 
within a sphere of radius 2cr around the nucleation centre. 

Across the traces of fig. l(a-d) we observe a gradual, 
tandem appearance of densification and ordering. The 
same pattern is evident in the additional new datasets 
as in the reanalysis of the old data. We attempt to make 
a local measurement of translational order distinct to 
orientational order by showing the standard deviation 
of the distance to nearest neighbours SD{r 2 n) over 


the first two shells for each particle 1(d). We make 
this direct-space measurement of translational order 
because the usual definition of translational order as 
relating to the presence of higher order Bragg peaks is 
inherently non-local and therefore difficult to apply to 
a small region of emerging solid phase, as well as being 
unnecessary in the case that direct-space information 
is available. Using this information we can see that 
translational and rotational ordering are synchronised, 
thus eliminating a candidate two-step nucleation mech¬ 
anism. By measuring the number of particles in the 
first neighbour shell via the algorithm of van Meel et 
al. [31] we show that the increase of density and order 
arises by an increase towards 12 in the average number 
of immediate neighbours 1(e), rather than by isotropic 
compaction or by expulsion of a thirteenth neighbour. 

All observables apart from the size of the largest 
discrete solid-like cluster start to smoothly deviate from 
their liquid values at t ~ —a^/6Di^ where Di is the long 
time self diffusion constant. The cluster size lifts off 
more sharply and a little later, because it is defined via a 
discrete thresholding of the q6 values. In refs. [6, 32] this 
threshold was set close to values expected for the bulk 
crystal. We suggest that the emphasis on densification 
as the leading process in ref. [6] arose mainly as result 
of this decision to treat crystalline order as a binary 
quantity, while it is now clear that local ordering varies 
continuously on the pathway between the fluid and 
crystalline phases. 

A neighbourhood of twelve particles is compatible 
with the (non-spacefilling) icosahedral symmetry as 
well as with fee and hep. There has been considerable 
discussion, especially in relation to glass formation, of 
the idea that this can lead to a multistep or arrested nu¬ 
cleation process, for a review see [33]. A small decrease 
in w6 (together with a large increase in q6) is associated 
with hexagonal ordering, however a larger decrease 
in w6 (together with a small decrease in g6) signifies 
icosahedra. An example threshold from the literature 
for icosahedral order is w6 < —0.023 [34]. By computing 
the w6 parameter as well as the q6 measurements we 
show that there is no signal of icosahedral order as 
an aspect of the nucleation process at the non-glassy 
overcompressions studied here (fig. 2). 

To further examine the time-correlation of the density 
and order changes, we plot single traces (the first three 
MC runs made at = 1.027, fig. 3). In the same 
way as for the average pathway, for individual runs the 
density and order increase simultaneously. Fluctuations 
evident in the individual traces, whether up or down, 
are correlated between the two observables with no dis¬ 
cernible time lag. From this we confirm the suggestion 
of fig. 1 that the increases of density and order are 
consequences of the exact same particle displacements, 
in contrast to the constant-density ordering which has 
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FIG. 2. The w6 spherical harmonic order parameter, used to 
identify icosahedral structure, does not show sufficient change 
(i.e. w6 < —0.023) to indicate that icosahedral order plays a 
role in nucleation at number densities (j) = [1.012.. 1.03]. 


been hypothesized [9]. 



FIG. 3. Time series of average qeqeii^i) and number density 
(inverse Voronoi volume, 0) for three MG runs at = 1.027, 
with averages collected in a sphere of radius 2a centred on 
the nucleation event at t = 0. Order and density increase 
simultaneously in each case, and fluctuations in the order and 
density traces are correlated with no time lag, see in particlar 
the upward-trending parts of each trace. The dashed line 
shows q6 = 0.28. 

In order to estimate the length and entropy scales 
associated with the initiating fluctuation, we plot the q6 
and number density averages leading up to nucleation 
for the set of MC simulations at = 1.027, with 
averages collected over spherical shells of increasing size 
around the point of initiation of the nucleus, defined as 
in fig. 1 (fig. 4, lines). For comparison we also map the 
log probability of observing a value higher than a given 
q6 or (j) in the bulk liquid, i.e. the cost in free energy 
associated with creating an order parameter fluctuation 
of this magnitude and radius. The trace at first final 
appearance of the solid-like nucleus (fig. 4, line with 
open circles) is associated with a free energy cost of 


approximately 3 to Gradual and simultaneous 

densification and ordering is indicated, with ordering 
initiated non-locally in a (potentially irregular) region 
of approximately 3 — 5cr radius. That the initiating 
density fluctuation does not become negative over the 
range plotted indicates that the inward mass flux is 
compensated at fairly long range. A version of fig. 4 
confirming similar behaviour at = 1.012 is appended 
as supplementary data (fig. S2 [30]). 



radius [cr] 


FIG. 4. (a) q6 and (b) inverse Voronoi volume, 0, averaged 
over spherical shells of increasing radii. The colour maps show 
the log probability to find a given order parameter value (or 
greater) when averaging over random shells of the same ra¬ 
dius in the bulk fluid. Traces indicate orientational averages 
over shells centred on the time-zero of the initiation of the 
nucleation event {first final solidification of a particle). The 
trace shading from light to dark indicates a time range of 
a^tjGDi = —1..1, with the t = 0 trace highlighted using cir¬ 
cles. The (pa^ = 1.027 MG dataset is shown here. Dashed 
lines indicate quiescent averages. 


III. CONCLUSION 

We have analysed the early stages of crystal nucleation 
in hard spheres. We observe that densification does 
not occur prior to bond orientational ordering, and also 
that bond orientational ordering does not occur prior 
to densification. Hard sphere nucleation starts with a 
process that includes densification and ordering. As far 
as we are able to measure orientational and translational 
ordering as distinct phenomena, we find that they also 
occur simultaneously. 

The early-stage hexagonally ordered fluctuations 
which we discuss here are identifiable with the ‘dense 
amorphous regions’ or ‘low symmetry clusters’ used by 
Schilling et al. [6] to argue the density-first case in 2010, 
and also appear at least similar to the ‘precursor struc¬ 
tures’ suggested to form without density enhancements 
by the Tanaka group [9, 13, 35] and used to argue an 
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order-first position. The explanation which we can offer 
for the historical divergence between density-first and 
order-first opinions is that probes of differing sensitivity 
to the respective phenomena have previously been 
employed. By showing that both variables lift from their 
quiescent distributions not only together on average, 
but with instantaneously correlated fluctuations in in¬ 
dividual traces, we hope that we have settled this dispute. 

We assert that nucleation of hard sphere crystals at 
low to medium overcompression begins with a collective 
fluctuation of radius approximately 4a, simultaneously 
manifested in the density and in the positional order. 
The radius of 4cr corrresponds to order 100 particles 
(depending on the shape of the nascent dense ordered 
region), indicating fluctuations of a highly collective 
nature. The increase of density at the nucleation site 
is supplied by long-range collective mass transport, 
with the amplitude of both density and order falling off 
gradually over the characteristic distance. 

We accept the current speculation that it is possible 
to advance a description of the hard sphere nucleation 
process as two-step, but only in the weak sense that it 
begins with an order/density fluctuation which initially 


strengthens more than it grows (before starting to grow 
more than it strengthens). At least when away from the 
glassy regime, the changeover between these two ‘steps’ 
is smooth and there is no sign of any qualitative changes 
in the type of order manifested. 

The most agreement that we can find with and across 
the recent literature is that nucleation is initiated with 
a diffuse entity of finite spatial extent: this model 
contrasts with the naiVe image of a sharply defined 
region of daughter phase expanding outwards from a 
point, and is a good candidate to eventually supply 
control and understanding of nucleation phenomena by 
matching the lengthscale and structure of experimental 
probes with those of the initial fluctuation. 
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